/*
 * LCells.h
 *
 *  Created on: Dec 18, 2014
 *      Author: marchi
 */

#ifndef LINKEDCELLS_SRC_LCELLS_H_
#define LINKEDCELLS_SRC_LCELLS_H_

#include <vector>
#include "Ftypedefs.h"
#include "Atoms.h"
#include "Metric.h"
#include "MyUtilClass.h"


using namespace std;
using namespace MATRIX;
using namespace DVECT;

class LCells {
	typedef vector<int> vectint;
	static double Rcut;
	static double Rcut_in;
	static vector<int> nc;
	vector<vectint>  indx;
	vector<vectint> Cellp;
	vector<vector<vector<vectint> > > Chainp;
	Matrix co,oc;
	vector<Dvect> x;
	int nr;
	vector<vector<int> > nnl;
	double Dist_ijk(int, int,int);
	void Init(Matrix & co0, vector<Dvect> & y){
		x=y;
		co=co0;
		oc=co.Inversion();
		nr=x.size();
		if(nc[XX] < 0){
			nc[XX]=static_cast<int>(co[XX][XX]/Rcut);
			nc[YY]=static_cast<int>(co[YY][YY]/Rcut);
			nc[ZZ]=static_cast<int>(co[ZZ][ZZ]/Rcut);
		}
		if(Chainp.size()) Chainp.clear();
		Chainp=vector<vector<vector<vectint> > >(nc[XX]);
		for(size_t m=0;m<nc[XX];m++){
			Chainp[m]=vector<vector<vectint> >(nc[YY]);
			for(size_t n=0;n<nc[YY];n++)
				Chainp[m][n]=vector<vectint>(nc[ZZ]);
		}
		if(nnl.size()) nnl.clear();
		nnl=vector<vectint>(nr);
		if(Rcut < Rcut_in) Rcut=Rcut_in;
	}
public:
	LCells();
	LCells(Matrix & co0, vector<Dvect> & y): nr(0) {
		this->Init(co0,y);
		this->Index();
	}
	static void Setup(double rcut, int nx, int ny, int nz){
		Rcut=rcut;
		int n[]={nx,ny,nz};
		nc=vector<int>(n,n+DIM);
	}
	static void Setup(double rcut){
		Rcut=rcut;
	}
	void operator()(Matrix & co0, vector<Dvect> y){
		this->Init(co0,y);
	}
	vector<vector<int> > & getNeigh(){
		try{
			if(!nnl.size()) throw " Neighbour list not yet computed !";
		}
		catch(const char * s){
			cout << s << endl;
			exit(1);
		}
		return nnl;
	}
	void Index();
	vector<vector<int> > & List();
	virtual ~LCells();
};

#endif /* LINKEDCELLS_SRC_LCELLS_H_ */
